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Complex eigenvalues, resonances, play an important role in large variety of fields 
in physics and chemistry. For example, in cold molecular collision experiments and 
electron scattering experiments, autoionizing and pre-dissociative metastable reso¬ 
nances are generated. However, the computation of complex resonance eigenvalues is 
difficult, since it requires severe modifications of standard electronic structure codes 
and methods. Here we show how resonance eigenvalues, positions and widths, can 
be calculated using the standard, widely used, electronic-structure packages. Our 
method enables the calculations of the complex resonance eigenvalues by using an¬ 
alytical continuation procedures (such as Pade). The key point in our approach is 
the existence of narrow analytical passages from the real axis to the complex energy 
plane. In fact, the existence of these analytical passages relies on using finite ba¬ 
sis sets. These passages become narrower as the basis set becomes more complete, 
whereas in the exact limit, these passages to the complex plane are closed. 

As illustrative numerical examples we calculated the autoionization resonances of 
helium, hydrogen anion and hydrogen molecule. We show that our results are in an 
excellent agreement with the results obtained by other theoretical methods and with 
available experimental results. 
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I. INTRODUCTION 


A. Motivation 


One of the biggest challenges of electronic structure calculations is to take autoionization 
into consideration. Autoionization is a process in which an electronic metastable state 
decays through spontaneous emission of an electron. Such a state has a finite lifetime, and 
is known as a resonance. The possibility that in a specific geometrical structure the molecule 
can be autoionized implies that the electronic and the nuclei coordinates are coupled. In 
such a case, the Born-Oppenhiemer (BO) approximation breaks down completely, making 
the electronic structure computation much harder. A possible solution to the problem 
is to impose outgoing boundary conditions on the eigenfunctions of the time independent 
electronic Hamiltonian within the framework of the BO approximation. In this way, complex 
potential energy surfaces (CPESs) are obtained. The real and the imaginary parts of the 
CPES provide, respectively, the energy (position) and the autoionization decay rate (width 
or inverse lifetime) of a molecule as function of its geometry. CPESs can be obtained by using 
complex basis functions jl|, analytical continuation of the Hamiltonian’s matrix elements 
Q or one of the complex scaling transformations, such as the uniform Q, exterior JJ] or 
smooth exterior scaling [5]. Alternatively, CPESs can be obtained by introducing a complex 
absorbing potential (CAP) jf], or a reflection-free CAP (RF-CAP) to the Hamiltonian £]• 
Essentially, the complex electronic eigenvalues obtained within the BO picture serve as 
potentials of the nuclear time-independent Schrodinger equation. In other words, by using 
the complex electronic eigenvalues there is no need to go beyond the BO approximation, 
since they introduce couplings between the nuclear and the electronic coordinates in a simple 
way. Still, although this idea was presented many years ago Jf|, so far it did not become 
a useful approach in the dynamics study of polyatomic molecules with many electrons. 
Unfortunately, the reason for this lies in the fact that most of the commercially used codes 
do not support the above modifications, rendering the computation of complex potential 
energy surfaces an unconventional task. However, the need for reliable CPESs is a must 
in large variety of fields in chemistry and molecular physics, in which quantum mechanical 
dynamics of molecules is of interest. A good example for this acute need is the most recent 
cold chemistry 0 and electron scattering experiments jllj]. Recently there were serious 
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efforts to develop codes for calculating CPESs (g, 12 16], Yet, it is most desirable to have a 


simple approach that uses standard electronic structure codes in order to calculate CPESs, 
since these codes are highly optimized and very efficient. 

Here, we propose a simple method of calculating CPESs, which utilizes standard elec¬ 
tronic structure codes without modifying them. This approach is based on an analytical 
continuation of results obtained in the real space into the complex plane. There are many 
different approaches to carry out analytical continuation from the real to the complex plane. 
Several of which are briefly described in the next section. In our approach a single real 
eigenvalue obtained from standard stabilization calculations is analytically continued into 
the complex plane. This idea is not new and it faced criticism in the past since the whole 


eigenvalue plot is not an analytical function of the scaling parameter 


& 


The transition 


from a stabilization plot into the complex plane and the search for a stationary resonance 
state goes through a singularity point, known as a branch point (BP) [18]. However, the 
way we implement this idea avoids these problems by not using the whole stabilizations plot 
but only an analytical part of it. 


B. Background 


As mentioned above, complex energy surfaces are the electronic eigenvalues obtained 
within the BO picture. Their real and imaginary parts correspond respectively to the energy 
position and width. These eigenvalues represent metastable states with a finite lifetime, res¬ 
onances. According to the Balslev and Combes theorem, atomic autoionization resonances 
becomes square integrable functions by applying uniform complex scaling transformation 
of the form r —> rrj, where r is the electronic coordinates and // is the scaling parameter 
defined as rj = aexp(i9) (and a and 9 are real) [ 3 !, 19, 20]. In order to implement Balslev 
and Combes theorem in a molecular system, Simon proposed the use of an exterior scaling 
transformation that avoids the singularities in the BO molecular Hamiltonian j^J]. In both 
atomic and molecular transformations the Hamiltonian’s spectrum is changing: while the 
bound states are unaffected and are characterized by real eigenvalues, the continuum states 
are rotated into the complex plane by an angle of 29, i.e., these complex eigenvalues strongly 
dependent on 9 [5]. In addition, the autoionization resonances, characterized by complex 
eigenvalues, appear in the spectrum. The first time they are exposed is at a critical value of 
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the rotation angle 6 = 9 bp- At this critical value a branch point is obtained in the spectrum 


of the complex scaled Hamiltonian (see Fig. 3 in Ref. 


2d, 


21|). However, as 6 increases the 


autoionization resonances stabilize and show low dependency on 0, meaning that stationary 
points in the complex energy plane are obtained. 

There are a few early approaches for calculating the resonances position and width from 
standard hermitian electronic structure calculations. For example, resonances can be ap¬ 
proximated from the density of states in the continuum, they can be calculated from the 
asymptotes of the continuum eigenfunctions or from the phase shifts of the eigenfunctions 


(see chapter 3 in Ref. [sj]). Alternatively, the resonances can be obtained from stabilization 
calculations where the eigenvalues are computed for an increased number of basis functions 


22 ] or when a finite given basis functions are scaled by a real factor 


23j, 


Here propose to analytically dilate the real energies into the complex plane via the Pade 
approximant. At its basis this approach is also not new: In 1981 Simons used the stabi¬ 
lization calculations and suggested to carry out a unitary transformation from the adiabatic 


energy levels to the diabatic presentation 


25]. In the diabatic representation the electronic 


energy levels are coupled, unlike in the adiabatic picture, where the electronic energies are 
represented as non-interacting states. Therefore, in the adiabatic picture avoided cross¬ 
ings are obtained, while in the diabatic picture crossings are exposed in the complex plane. 
These crossings represent the branch points, and their exposer facilitates locating the nearby 
stationary points (the resonance energy). 


An even simpler approach was introduced by Thompson and Truhlar 


26] and later mod¬ 


ified by Isaacson and Truhlar (27|. Under this approach the resonance complex energy is 
obtained by analytically continuing a single eigenvalue into the complex plane, to this end 
the whole stabilization curve is used. However, as McCurdy and McNutt 17] pointed out, 
analytic continuation of a single-root fails due to the existence of non-analytic regions in 
the eigenvalue plot. Hence, McCurdy and McNutt suggested a multi-eigenvalue method. 
They carried out analytical continuation of the characteristic polynomials of the Hamilto¬ 
nian matrix using at least two eigenvalues. Doing so, they avoided the non-analyticity of the 
corresponding eigenvalues and located stationary points. In fact, McCurdy with McNutt 
opened a new research direction for calculating resonances from stabilization graphs. This 
multi-eigenvalue method is indeed pursued until today j^S, 3] and it relies on the correct 
description of the avoided crossings. 
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Nevertheless, in this paper we show that under easy to fulfill conditions, resonances can 


be calculated from a single stabilization root, similar to Thompson and Truhlar 


26J. Yet, 


unlike Thompson and Truhlar, we suggest to analytically dilate the eigenvalue as function 
of the stabilization parameter while avoiding the branch point areas, i.e., the non-analytic 
structure of the complex energy. This can only be achieved when carefully choosing the 
area in the stabilization plot that will be dilated, meaning, excluding the avoided crossings. 
In other words, we are only using the stable analytical region in the stabilization graph. 
We are not interested in describing branch points. Furthermore, we illustrate later that 
we can always remain in an analytic area, and eventually converge to a stationary point in 
the complex plane. The converged energy is in an excellent agreement with the eigenvalue 
obtained by an explicit complex scaling calculations. 

The proposed analytical continuation scheme resembles Moiseyev’s and Corcoran’s sug¬ 
gestion to treat molecular resonances within the BO picture by evaluating the molecular 


Hamiltonian matrix elements 0]. These matrix elements are analytical functions of the 
scaling factor, even though the operator is not, because the contour of their integration in 


the complex plane can be chosen to be such that avoids the singular points tl, B, 0-|32l. 
Therefore, they can easily be analytically continued to the complex plane [1, 2 }]. In a similar 
way, we demonstrate here the existence of an analytical path from the stabilization graph 
towards a complex stationary point. This path bypasses any branch point, thus avoiding 
any singularity. 


II. THEORETICAL SCHEME 


A. Resonances from square integrable basis functions: complete vs. finite bases 


Upon complex scaling, autoionization resonances become square integrable functions, 
therefore are exposed as eigenfunctions in the Hilbert space. When the complex scaling 
parameter is defined as r/ = aexp(id), the continuum states are strongly dependent on 0, 


however, in the exact limit, they are invariant to a 
set, the continuum states stay fixed as a is varied 


• T 


rat is, when using a complete basis 


12-i 


Contrary, in actual calculations, when incomplete and finite basis sets are employed, the 
continuum spectrum does vary with respect to a. An example for this behavior is seen in 
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Fig.ffl In this figure the relevant energy levels of the helium atom and the hydrogen molecule 
are plotted as a function of a. In both cases a finite number of Gaussians were used as a basis 
set in the calculations (further computational details are given below) and it is clear that 
the states depend on a. This dependency can be utilized in order to calculate the resonance 
position and width. To do so one must first identify the resonance footprints in the real 
energy spectrum. These footprints lie in the stabilization plot. In Fig. Q] these footprints 
are clear: while continuum energies change strongly with respect to a, one specific energy 
seems to stay fixed and emerges out of the plot. This energy corresponds to the resonance 
position which is not highly affected by the variation of the scaling factor. The reason 
for this stabilization is clear: the resonance wave-function is localized in space. Therefore, 
relatively small expansion or compression of the basis functions does not have much effect 
on it. In contrast, the other states are associated with delocalized functions and therefore 
are strongly affected by the scaling factor. 

Although the resonance position emerges out of the stabilization plot, its width does not. 
One way to obtain its width is to move into the complex plane. There, the resonance is square 
integrable, and appears as a complex eigenvalue of the Hamiltonian. The imaginary part 
of this eigenvalue corresponds to the resonance width. One easy method to move into the 
complex plane is by analytical continuation from the stabilization plot jl2, 2^3- However, 
this continuation should be done with caution since an avoided crossing in the stabilization 
plot is associated with singularity, a branch point in the complex plane. This singularity is 
caused due to the fact that as the resonance becomes square integrable, it separates from 
the continuum. At this critical point, the resonance and the rotating continuum coalesce in 
the complex plane. This causes the spectrum to be deficient and thus non-analytic. 

The existence of an avoided crossing in the stabilization graph raises difficulties, and 


analytically continuing a single eigenvalue can fail to accurately 
the complex plane 17]. For this reason, McCurdy and McNutt 


coworkers 


28,|29] 


ind a stationary point in 
17| and later Jordan and 


mm used at least two eigenvalues to obtain a truncated characteristic polyno¬ 


mial of the Hamiltonian. Then, they performed analytical continuation of the polynomial 
coefficients instead of the eigenvalue. Yet, since finite basis sets are always employed, there 
is a way to analytical continue a single eigenvalue into the complex plane. The “trick” is 
to use only a small part of the stabilization plot, a relatively stable region, instead of the 
whole plot. This region is an analytical function of a and it avoids the avoided crossings, 
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therefore, it is a perfect starting point for the desired analytical continuations. 

It is important to stress that this “trick” is only applicable in a finite basis set. Under 
an infinite basis set it is impossible to carry out analytical continuation of an eigenvalue 
into the complex plane. As the basis set approaches completeness the number of avoided 
crossings increases, and thus the singularity area in the complex plane increases too. That 

is, there is no way to analytically dilate a single eigenvalue and avoid the branch point. 

In the next section, a comprehensive discussion and illustration are presented explaining 
why when finite basis sets are used, it is possible to bypass the singularity and find an 
analytical route towards the desired stationary point. Moreover, it is shown that when finite 
basis sets are used, it is actually very difficult to find a branch point let alone pass through 

it. 
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(a) (b) 

FIG. 1: [Color online] Energy stabilization plots for (a) the helium 2s 2 resonance state (b) the 
hydrogen molecular lcr 2 resonance state (R = 1.4 a.u.). The stabilization is obtained by varying 
the real scaling parameter, a. The red and green areas in (a) represent the two regions that were 
used as a starting point for the two analytical continuations. For H 2 only one stable region, marked 
in red in (b), was used for the analytical continuation. 
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FIG. 2: [Color online] Complex eigenvalues of a non-hermitian helium Hamiltonian matrix obtained 
by the use of uniform complex scaling, r] = <aexp(z0) where 9 is varied and a is held fixed. Three 
6 -trajectories are presented with different a values. The app trajectory is presented in black and 
includes also the branch point (pink). The green arrows represent the direction of the ^-trajectories 
from zero by 0.001 radians increments. Note, that the rotational angle 9bp •> in which the branch 
point is obtained, is much smaller than the value of 6 at the stationary point. In these uniform 
complex scaling calculations a finite Gaussians basis set is used. 
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B. A finite basis set as a means to bypass the branch point 


As previously discussed, the resonance position and width can be obtained by moving 
into the complex plane using the complex scaling parameter r/ = aexp(i9) . In the complex 
plane, the resonances are recognized as complex energies that are not affected by a small 
change in a and 0, i.e. stationary points. Graphically, we identify them as cusps in the 6 
or «-trajectories plots. Such stationary points are located close to the branch point, but 
are always deeper into the complex plane (see Fig. 2 in Ref. [l2] and Fig. 2 in Ref. 18)h 
Another good example for this is given in Fig. [2], where Gtrajectories calculations for the 
2s 2 resonance state of helium are presented. In these trajectories, a is held fixed and 
the rotational angle, 0, is increased from zero (the hermitian case). As indicated in the 
black ^-trajectory, the rotational angle 6 = 6bp, in which the branch point is obtained, 
is much smaller than the value of 0 at the stationary point. That is, in order to locate 
a stationary point by varying 6 from the real axis, it is necessary to go through a branch 
point. Theoretically, this feature is a problem when one tries to find the stationary point 
through analytical continuation, however as clearly seen in Fig. El in practice it is not. In 
this figure, two additional trajectories are plotted. These blue and red trajectories reach the 
stationary point smoothly without passing through a branch point. In other words, contrary 
to the complete basis set case, in which the branch point is a-independent, when a finite 
basis set is used, the branch point is obtained for a very specific value of a = asp (as shown 
in Fig. [ 2 ). As a matter of fact, it is very hard to calculate a branch point and there are 
specific methods developed in order to locate it jl8, Tij]. Here for example, the branch point 


parameters, 6bp and asp, were found using such a method [33]. However, these parameters 
are never fully exact but rather a mere estimate to the real values. To sum it up, for finite 
basis sets, the fact that the branch point is obtained only for very specific values of a and 
0 is an advantage. 

In light of the above, the resonance stationary point can be found via analytical con¬ 
tinuation of the real eigenvalue from the stabilization plot into the complex plane. This 
is examined below by comparing atomic and molecular resonance positions and widths ob¬ 
tained by the proposed analytical continuation with the results obtained by explicit uniform 
complex scaling methods. 

Below a mathematical demonstration is given for the claim that in a finite basis set 
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analytical continuation to the complex plane is feasible. It is based on the fact that a is 
always different than the exact a B p. 

Sufficiently close to the branch point where the stationary point lies (a = asp), the 


complex eigenvalues associated with t 


leading terms in the Puiseux series 


B 


ie resonance and continuum states are given by the 


m- 


E± = E bp ± b^/rj^rjBP 

= Epp ± by/aspe — ~a B p^^ p 

= E bp ± b^( - 1 ) , (1) 

\C*BP J 

Here, coalescence of two solutions can only occur when - 

A(0) = ^L e ie ~ ie BP _ i o, (2) 

OLbp 

Therefore, it is clear that a coalescence is an unlikely scenario, particularly when 6 = 0 B p- 
That is, since asp / a .bp the eigenvalue associated with the stationary point is an analytical 
function of 9. 
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C. Analytical continuation via the Pade approximant 
— from real eigenvalues to the complex energy plane 


In this paper, we analytically dilate a hermitian eigenvalue to the complex plane by the 
Pade approximant. More specifically, we generate an analytical approximation to E(//), 
where rj = a exp (id), by the Schlessinger point method This method requires a set 

of M input variables and their corresponding values. In our case, we select rj as the input 
variable, or more precisely rj = a (where 0 = 0), and E(a) as the corresponding value. The 
input data, E(a ), is taken from the relatively stable region in the stabilization plot which 
excludes the avoided crossing. For example, in order to calculate the 2s 2 resonance state of 
helium, we used the red or green regions in the stabilization plot shown in Fig. [TJ a). That 
is, each region (red or green) is used as a starting point for a different calculation. The 
Schlessinger truncated continued fraction has the form - 



( 3 ) 


where the Zi coefficients are chosen such that 


C M (oti ) = E(on), * = 1,2,..., M. 


( 4 ) 


Once the Zi coefficients are determined, we perform an analytical continuation into the 
complex plane by evaluating Cm(v) (where // is complex, i.e. 9 ^ 0). Convergence of the 
extrapolated function, Cm(v ), with respect to M is routinely checked, and the difference 
between Cm(v) an d Cm-i(v) is reported as the Pade error. 

The procedure we propose is summarized as follows: 

1) Calculate a stabilization plot using a standard electronic-structure approach as in Fig. Q3 

2) Use the points in a stabilized region as input for the Schlessinger point method. 
Remember to exclude the avoided crossings like the colored regions in Fig. [lj 

3) Fit the selected points to a Pade polynomial by making sure that this polynomial 
reproduce the original values for the input variables. 

4) Preform analytical continuation in the real plane (i.e., with 9 = 0) to make sure an 
appropriate set of points was chosen. A too dense set of input points will damage the 
performance of the Pade approximant, and will end up in a linear extrapolated plot that is 
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almost independent of a. 

5) Preform analytical continuation into the complex plane: generate a ^-trajectory with a 
fixed a. The fixed a should be taken from its values in the stable region of the stabilization 
plot. 

6) Look for an optimal 6 value in which a cusp is seen in the 0-trajectory. 

7) Calculate an a-trajectory using the optimal 6 found in step 6. Look for the optimal a in 
which a cusp is seen in the cr-trajectory. 

8) Repeat steps 5-7 until the optimal a and 6 are converged. Upon convergence a clear 
cusp is obtained in both a and ^-trajectories. Both cusps will touch each other (see Fig. EO). 

lese cusps are associated with a stationary point that satisfy £-E(rf) —>■ 0 and E(r) ) —> 0 


a. 


9) Check the Pade errors of the stationary point and make sure it is reasonable. 

10) It is recommended to repeat the process with a slightly different set of input points to 
examine the stability of the stationary point. 


This procedure is very easy and quick. The only time and computational consuming 
step is step 1. That is, this procedure is much simpler than actual electronic-structure 
calculations inside the complex plane. 


III. COMPUTATIONAL DETAILS 


In this work we study resonance states of the helium atom, the H _ anion and the H 2 
molecule with the leading 2s 2 , 2s 2 and lcr 2 configurations, respectively. In all cases we reach 
the full configuration interaction (FCI) limit. For helium we use our own Cl code for two- 


electron atoms and ions including complex scaling transformation 


35]. For H and H 2 we 


use the standard equation-of-motion couple-cluster with singles and doubles (EOM- CCSD) 


for excitation energies implemented in the Q-Chem package 


36]. The ground states, which 


play the role of the reference states, are Is 2 for helium and H , and la 2 for H 


12 - 


The basis sets employed comprise of primitive Gaussian functions since they can be scaled 
trivially. Uniform complex scaling is performed by rotating the Gaussians basis functions 
by r —» rry, where r is the Gaussian coordinate and rj is the scaling parameter defined as 
rj = a exp(iO) (a and 9 are real) jd, fd], 2tJ. Stabilizations are obtained by scaling the 
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Gaussians when 0 is set to zero. 

For helium we use a 19s 15p 10 d 8/ (19/15/10/8) series of exponentially tempered basis 


functions 


35], The basis set is optimized for calculating the helium states with the highest 


principal number n —2 and the highest partial angular momentum 1=3 (Is 2 , ls2s, 2/ 2 ) 
for about 10 -5 Hartree precision. 

For H - we use a 12s lOp 5 d 3/ 2 g (12/10/5/3/2) set based on the Bentley and Chipman 
basis, which include s- and p-type functions, and was used to calculate the H _ resonance 


before 


0 - 


The contracted Is orbital is transformed into three primitive functions. In 


addition, the basis is augmented with the d-, f- and (/-type functions of the aug-pcJ-4 basis 
set j^l . Convergence of the resonance energy with respect to the number of diffuse basis 
functions is ensure by further augmenting this basis with two s- and two p-type diffuse 
functions, which results in a 14/12/5/3/2 basis set. 

For H 2 we report results obtained by using a 12s lOp 3d 2/ (12/10/3/2) basis, which is 


[33]. 


based on the 5ZP basis set 139(1. In this basis we replace the most diffuse s-function with 
six diffuse functions, and the most diffuse p-function with seven diffuse functions, where an 
even-tempered spacing of 2.0 was employed for their construction. In addition, the most 
tight d- and /-type functions where replaced with diffuse functions. This basis is referred to 
as Basisl. Basis set convergence is ensured by comparing this basis set results with two other 
bases. First, to make sure we obtain convergence with respect to the number of diffuse basis 
functions, we compare the 12/10/3/2 results with the ones obtained with a similar basis in 
which the most diffuse s and p functions are removed (a 11/9/3/2 basis set). We observe 
convergence within 10 -8 Hartree for the resonance energy between the two bases. Second, 
we examined a denser 12/9/3/2 basis set, in which the three most diffuse s -functions of the 
original 5ZP basis set are replaced with eight diffuse functions (with even-tempered spacing 
of 2.0), everything else is similar to the 11/9/3/2 basis set. This basis is referred to as 
Basis2. Below, we report results obtained by Basisl and Basis2. 

All the bases exponents are presented in the supplementary information. 
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IV. NUMERICAL APPLICATIONS 


A. Helium autoionization resonance — a comprehensive test case 


Since helium is a well-studied and simple system, 


35( the helium 2s 2 state can serve 


as 


a comprehensive test case, for which we can study autoionization processes and explore the 
capabilities of new computational schemes. Furthermore, helium is a two electron system, 
hence it is easy to calculate this resonance position and width using FCI and uniform complex 
scaling (UCS). Doing so, no approximation is done on the electronic structure and a pure 
comparison to a new method can be performed. Therefore, in our study this calculation 
was used as a reference point. Such a comparison for helium 2s 2 resonance can be seen 
in Fig. El In this figure, optimal cusps from both UCS calculations and our analytical 
continuation calculation are presented. For the analytical continuation calculations the 
cusps were obtained through the iterative procedure described in Section III Cl A similar 
iterative scheme, in which FCI calculations were performed in each step, was carried out 
for the UCS calculation. A remarkable agreement between the UCS calculations and our 
analytical continuation is clearly observed. 

In Fig. El(b) two cusps are displayed. One cusp was obtained by using region 2 in Fig Ufa) 
as the starting point for the Pade analytical continuation. The converged scaling parameters 
of the Pade cusp were then used as starting points for the UCS iterative calculations. Doing 
so, a UCS cusp was also obtained. The two cusps were in an excellent agreement, both 
in their scaling parameters and in their position and width. In other words, the Pade 
approximant was able to produce a very accurate a- and ^-trajectories based on the real 
stabilization plot. 

In Fig. Etc) two cusps are displayed for the Pade analytical continuation. This time, the 
cusps were obtained using region 1 in FigEta) as the starting point for the Pade approximant. 
Each cusp represented another resonance stationary energy, where the two energies differ 
by 2 x 10 -4 Hartree for the imaginary part, and by 4 x 10 -5 Hartree for the real part. 
Similarly, in Fig Etd) two cusps are also displayed for the UCS calculations. Again, each 
cusp represented another resonance stationary energy where these two energies are within 
less than 10 -4 Hartree difference for both real and imaginary parts. It is important to 
note that the UCS cusps were obtained as described above. That is, the converged scaling 
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parameters of the Pade cusps were used as starting points for the UCS iterative calculations. 

A good summary of Fig. HR-d) is displayed in Fig. EH a). This figure demonstrates the 
proximity of the Pade results to the UCS ones. Each stationary point obtained by the UCS 
calculations has an analogues point obtained by the Pade approximant, where in all cases 
there is a good agreement between them. In fact, the distance between each cusp couple 
was calculated and was found to be: 1.4 x 10 -4 , 1.1 x 10 -4 and 2.5 x 10 -4 Hartree for the 
black, red and green cusps, respectively. Recalling the fact that the Pade procedure is a fast 
and simple computational scheme, these results are very encouraging. 

In order to better understand how our analytical continuation scheme works, a thorough 
investigation was done using the UCS calculations. During this investigation, the complex 
eigenvalue associated with He (2s 2 ) resonance was calculated at 720 different scaling pa¬ 
rameters (9 and a). Fig. H represents the absolute value of this eigenvalue derivative at 
these points in both 3D and 2D contour plots. Fig. Ha) and (c) display this derivative with 
respect to 9, whereas Fig. HR) and (d) display this derivative with respect to a. Look¬ 
ing at the figures, the theory behind our analytical continuation becomes clear. First, it 
is obvious that for small 9s, there are regions in the derivatives (with respect to a and 9) 
that exhibit very large values. These values are associated with the branch points and the 
avoided crossings, which are forbidden areas for analytical continuations. Yet, it also quite 
clear that in spite of these forbidden areas, there are analytical paths that start on the real 
scaling parameter axis (i.e. 9= 0) and end up in stationary points in the complex plane. 
Moreover, there are even paths that lead from one resonance stationary solution to another. 
An example of such a path is the gray dashed line in Fig.H( c ) and (d). This path starts from 
region 2 in the stabilization plot of Fig. Ufa), and goes through each stationary point that 
was found in Fig. [HI That is, this investigation is a numerical illustration for the possibility 
to find analytical paths from areas in the standard stabilization calculations to the complex 
stationary points. 

Another interesting point that can be seen in Fig. His that the entrance to the complex 
plane from the real axis always starts as a very narrow path. In fact, these passages occur 
only since the calculation employs a finite basis set. For an infinite basis set the passages will 
close down, since as the basis set approaches completeness the number of avoided crossing 
increases, and the singularity area in the complex plane increases too. That is, there is no 
way to analytically dilate the single eigenvalue and avoid the branch points, the passes are 
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close. 

The existence of only narrow paths in a finite basis set framework, indicates that a careful 
search for these analytical paths in the complex parameter plane (a and 9) is required and 
still, the Pade approximant has no problem locating these paths. However, in the absence 
of such analytical paths the method presented here is not applicable and one may use the 
analytical continuation of the characteristic polynomials as described in Ref. fl7 ]. 

Fig. 0 is another illustration to the information in Fig. HI In this figure we show a 3D plot 
of the helium 2s 2 complex eigenvalue as 6 is varied, i.e., we show a-trajectories at different 
fixed 9s. In this figure we examine the complex eigenvalue itself, unlike in Fig. 0] where the 
derivative of the complex eigenvalue was examined. It is clear that cusps are obtained in 
the a-trajectory calculations as 9 becomes larger. However, it is also clear that there are 
paths connecting the different cusps and that these paths start from very specific regions. 
These regions are the stable parts of the a-trajectory at small 9s, where the eigenvalues 
are relatively close to one another. In other words, these paths go to the stable part of 
the stabilization plot. In Fig. [5] these paths are indicated in the warm colors (red-yellow), 
whereas in Fig. 0] such a path is marked by the gray dashed line. 
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FIG. 3: [Color online] Helium 2 s 2 resonance stationary points and the corresponding a- and 
0-trajectories, (a) Stationary points energies for the 2 s 2 resonance. Empty triangles represent 
stationary points obtained by UCS, while the full triangles are stationary points obtained by our 
analytical continuation. Each triangle correspond to the cusps in (b), (c) and (d). (b) a- and 
0-trajectories obtained from our analytical continuation scheme starting from region 2 in Fig. [l]( a) 
(green lines) and the corresponding UCS calculations (purple lines), (c) a- and 0-trajectories 
obtained by our analytical continuation scheme starting from region 1 in Fig. [l](a). (d) a- and 
0-trajectories obtained from the UCS calculations that correspond to (c). In (b) a clear cusp is 
seen, indicating one stationary point. In (c) and (d) two cusp are seen (one in black and the other 
in red) indicating two stationary points. Every stationary point is marked with a blue arrow. In 
(b) (c) and (d) the closed and open circles correspond to the a- and 0-trajectories, respectively. 
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e e 

(C) (d) 


FIG. 4: [Color online] 3D plot and the corresponding 2D contour plot of \d(E — iY)/d6\ (a, c) and 
| d(E — iT)/da\ (b, d) of the 2 s 2 helium resonance as function of a and 6. White arrows mark 
the complex stationary solutions, for which the complex derivative is minimal. These points are 
associated with the cusps in Fig. [3l Note the narrow analytical paths that go from the real axis to 
the complex plane. A gray dashed line demonstrates such a path. This path starts from a certain 
area in the real stabilization plot and goes to the complex plane through the three stationary 
points. 
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a trajectory 



FIG. 5: [Color online] A 3D plot of the helium 2 s 2 complex eigenvalue as 6 is varied. This plot 
shows a -trajectories at different fixed 0’s. The stationary solutions (cusps in the a trajectories) 
are clearly seen for larger values of 0, and it is clear that they can all be connected through certain 
paths. For rather small 6 a large dispersion in the complex eigenvalue is shown. However, the 
paths that connect the different cusps start from a relatively dense area marked by the warmer 
colors (red-yellow). These areas correspond to the stable part of the stabilization plot in Fig. 1(a) 
(the red and green region there). 
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B. Autoionization Feshbach resonance of the hydrogen atomic anion 

In order to further examine our Pade analytical continuation scheme we calculate the 
2s 2 Feshbach resonance of the hydrogen atomic anion. The results are presented in Table U 
where they are compared to experimental and other theoretical works. An excellent agree¬ 
ment with these estimations is observed, as our results are well within the experimental 
errors. Furthermore, the energy position in this work is lower by 2 x 10 -4 Hartree from other 


theoretical works, and the width is only 1.5x10 4 Hartree lower than Ho 


)’s (lO 


and Chen’s 


4JJ results, while it is 2x10 4 Hartree higher than the width obtained by Bravaya et al. 
12 | . Additional computational details and further results can be found in the supplementary 


information. 


TABLE I: The real part of the energy (ReE) and width (T) of the 2 s 2 Feshbach resonance of H 
in Hartree 


Ref. 

-ReE 

r x lo -3 

Experiment 

McGowan (1967)“ 

0.1485 

1.6 


±0.0004 

±0.2 

William (1976) fe 

0.1488 

1.65 


±0.0004 

±0.2 

Theory 

Ho (1981) c 

0.1487765 

1.731 


±0.000002 

±0.0008 

Chen (1997) d 

0.148782 

1.72 

Bravaya et al. (2013) e 

0.1488 

1.38 

Present work (Pade) 

0.14855 

1.56 


±0.00001 

±0.03 


a - Ref 


421] , & - Ref [43j, c - Ref (40|, d - Ref j4jJ, e - Ref 


n 
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TABLE II: The real part of the energy (ReE) and width (T) of the 1 Feshbach resonance of H 2 
(i?=1.4 a.u.) in Hartree, ReE is presented with respect to the ground state energy (-0.56994 
Hartree). 


Method 

ReE 

r 

Complex CI a 

0.4630 

0.0272 

CMCSCF a 


0.4638 

0.0270 

RF-CAP fe 


0.4615 

0.0227 

Basisl (Pade) c 

0.4619 

0.0230 




±0.0002 

±0.001 

Basis2(Pade) c 

0.4601 

0.0231 




±0.00001 

±0.0004 

a - Ref 

44 

]. b- Ref p. c- 

Present work. 


C. Autoionization Feshbach resonance of the hydrogen molecule 


Calculating molecular resonance is a complicated task, which poses a challenge on the 


regular complex scaling methods 


0 , 4, 


3CH32J. Therefore, it is important to test the Pade 


analytical continuation scheme presented here on a molecular system such as the hydrogen 
la 2 molecular resonance. This Feshbach resonance was calculated at an internuclear distance 


of R = 1.4 a.u. for which there are several calculations available for comparison. Note that 
in the molecular case, when the basis functions exponents are scaled by the real factor, a, 
the electronic coordinates of the Gaussian basis functions are shifted from their centers, A, . 
That is, every one of the Gaussian basis functions is scaled as, G(fi — Aj) —» G([r; — Aj]/a). 

In Fig. E an optimal cusp for hydrogen la 2 molecular resonance is shown using Basisl 
(see Section HTT1) . This cusp was obtained through analytical continuation of the stable region 
in the stabilization plot of this resonance (marked in red in Fig. [U(b)). 

Table HT1 presents our results using both Basisl and Basis2, where Basis2 is a denser basis 
than Basisl (see Section EU). We observed relatively minor effect when comparing these 
basis sets. The difference for the real part was in the order of 0.0018 Hartree whereas there 
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ReE, Hartree 


FIG. 6: [Color online] a-trajectory (black) and 0-trajectory (red) obtained from our analytical 
continuation scheme for the H 2 (R = 1.4 a.u .) 1 autoionizing resonance. In this figure an 
obvious cusp is seen (blue arrow), indicating a stationary point at 0=0.18125 and 0 :=0.658. The 
a-trajectory and 0-trajectory overlap at the cusp, as clearly shown in the inset. 


was no difference in the imaginary part. 

In Table [TT1 we compare our results with other theoretical works. Our results are in 


an excellent agreement with these works, particularly, with the RF-CAP jg]. A good agree¬ 
ment with the complex Cl and complex multi-configuration self consistent fields (CMCSCF) 


methods is also observed 


44], 


V. CONCLUDING REMARKS 

In this paper, we demonstrate that a single-eigenvalue curve taken from standard stabi¬ 
lization plots can be analytically dilate into the complex plane. On this basis, we suggest 
a simple procedure that accurately locates resonance positions and widths. Hence, it opens 
up a new way to calculate CEPS utilizing standard electronic-structure codes. This method 
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was successfully tested in calculating helium 2s 2 resonance state and was compared with 
explicit uniform complex scaling results yielding an excellent agreement. 

In addition, an in depth analysis was presented. In this context, numerical illustration 
clearly indicated the presence of analytical paths from the real axis to the complex plane 
when using finite basis sets. These paths start from the stable region in the standard 
stabilization plot. They continue through a very narrow passage surrounded by huge barriers 
that represent the avoided crossings and branch point areas. Finally, these paths end up in 
valleys of stationary points. 

Since these paths are analytical, they can be traced using the Pade approximant. In 


this work the Pade approximant was generated by the Schlessinger point method 


34J. The 


input points for this analytical continuation were taken from only the stable part in the 
stabilization plots. In this way, it is simple to infiltration through the analytical paths and 
find the resonance stationary points. It is important to note that this procedure fails when 


0 . 


input points are taken from the whole eigenvalue curve 

We showed that the success of the proposed approach is an outcome of our shortcomings: 
the mandatory use of finite basis sets in numerical calculations makes the existence of a 
branch point in the complex plane a rare occasion. In this case, analytical paths from the 
stabilization plot to the stationary point emerge. 

Finally, we implemented our approach on different chemical systems. In addition to 
helium, the resonance energy of the hydrogen 2s 2 atomic anion was calculated as well as 
the molecular Feshbach la 2 resonance of H 2 . The results for the hydrogen 2s 2 atomic anion 
were in an excellent agreement with experimental and other theoretical evaluations, whereas 
the results for the H 2 1 a 2 resonance were in an excellent agreement with other calculations 
performed explicitly in the complex plane. 

In spite of the above, we do not claim that our approach avoids the need to develop 
and use other methods. The reason for this is the fact that a region stable enough for 
a successful analytical continuation is not always guaranteed. It is recommended to first 
verify the conjecture of analyticity at certain geometries, particularly in the calculations of 
molecular autoionization resonances. This can be done by comparing the results with other 
methods. As the next step, one can use the method we proposed for calculating the entire 
complex potential energy surface. 

To sum it up, here we open a new possibility for calculating atomic and molecular autoion- 
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ization resonances in a very simple manner and substantially lower computational efforts. 
This method has the potential to bridge between standard electronic-structure methods and 
calculations of CPESs, which are highly important in studying the dynamics of a molecular 
system for example. 
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